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It is believed that first-order phase transitions at or around the GUT scale will produce high- 
frequency gravitational radiation. This radiation is a consequence of the collisions and coalescence of 
multiple bubbles during the transition. We employ high-resolution lattice simulations to numerically 
evolve a system of bubbles using only scalar fields, track the anisotropic stress during the process 
and evolve the metric perturbations associated with gravitational radiation. Although the radiation 
produced during the bubble collisions has previously been estimated, we find that the coalescence 
phase enhances this radiation even in the absence of a coupled fluid or turbulence. We comment 
on how these simulations scale and propose that the same enhancement should be found at the 
Electroweak scale; this modification should make direct detection of a first-order electroweak phase 
transition easier. 



Gravitational waves are a direct probe of the physics 
of the early Universe and the only direct probe of 
physics before recombination. The direct detection of the 
stochastic gravitational wave background will carry with 
it new insight to high energy physics and the nature of 
gravity. Cosmological processes that produce stochastic 
gravitational wave signals — inflation, cosmic string net- 
works, phase transitions — are of substantial interest since 
they represent an untapped source of information. 

Direct detection experiments such as the Laser In- 
terferometer Gravitational- wave Observatory (LIGO) 
are now underway and the hope of future missions, e.g. 
Q , will further the direct detection search; the time for a 
priori estimation of signatures will soon come to an end. 
Our current intention is to provide precision estimates of 
the gravitational wave spectrum to aid in the design of 
the next generation of observatories. 

The spectrum of gravitational waves from first-order 
phase transitions has been of interest for a number of 
decades. Although estimates of the gravitational radi- 
ation from cosmological processes date to the the early 
1980's [3j, the first rigorous attempts at estimating the 
power radiated from the collisions of bubbles came some- 
what later, first from vacuum bubbles 0, Q and then in 
the context of phase transitions 0, 0| • It was in this later 
work that the authors discovered that turbulence would 
additionally contribute to the stochastic gravitational 
wave spectrum. In recent times, interest in the gravita- 
tional radiation produced during phase transitions, with 
particular emphasis on the electroweak phase transition 
[8h12| , has led to good estimates and has rejuvenated the 
discussion. These estimates accurately address the power 
from detonations: highly energetic domain walls colliding 
at large velocities and a subsequent period of turbulence 
(the authors of [13| have additionally addressed the para- 
metric dependence of that terminal velocity to the results 
of the preceding work) . Analytic work has begun to ad- 
dress the added complexity of deflagrations. In this case 
thick walls (preceded by shock fronts) smear out the large 
field gradients, and exact descriptions of these collisions 



will need to be explored using fluid dynamics. We leave 
the topic of deflagrations to a future manuscript. 

Although much of this work has also addressed the 
added contribution of gravitational wave power from tur- 
bulence produced by the coupled fluid during the phase 
transition, we aim to present an additional effect. Even 
in the absence of proper fluid turbulence, a coalescent 
phase of the fields participating in the phase transition 
amplifies the gravitational wave signal and shift the peak 
frequency. Although this is discussed as possible for the 
electroweak phase transition, see e.g. [14j . we see this as 
a toy model where we concentrate on the generation of 
gravity waves from a scalar-only phase transition. 

Here, we discuss the possibility that the Universe un- 
derwent a first-order phase transition at or near the GUT 
scale, ~ 10 13 — 10 15 GeV, at time t*. Later this time will 
correspond to a program time, r = 0. We realize the 
phase transition with a scalar field, <f>, and a Coleman- 
Weinberg type associated potential, 



(1) 



where e parameterizes the height difference between the 
two minima, A is the dimensionless self-coupling of the 
fields and </>o sets the scale of the transition. We param- 
eterize the energy density, p*, of the whole Universe at 
the time of the transition by the energy scale, p,, such 
that 



(2) 



Just before the transition begins, the field is static, so 
the homogeneous energy density of the field is given by 
the energy of the metastable minimum, 



p<j,{U) = 2e\4>Q. 



(3) 



We presume that the field undergoing the phase transi- 
tion is some fraction, £, of the total content of the Uni- 
verse at that time, so the energy in the field is 



(4) 
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Most of the literature that studies phase transitions pa- 
rameterizes energy using the ratio of the energy of the 
false minimum compared to the overall thermal energy 
0, for this case, 



(5) 



The field obeys the usual Klein- Gordon equation in an 
expanding background, 



3iJ</>- 



V 2 </> dV(<f>) 



= 0, 
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the ratio of the bubble radius, r*, to the initial Hubble 
length, . We chose r^H^ 1 — 0.07 so that up to five 
bubbles can be easily nucleated on the lattice, where the 
initial length of each side is H^ 1 /2. We will always use 
a 3-dimensional lattice with finite-time differencing and 
we define program time r = t — t* . 

We begin each simulation by randomly placing N bub- 
bles on the grid. The average spacing between the bub- 
bles, L*, is related to the nucleation rate /3 and can be 
calculated by \2l\ 



.55396n- 1/3 fl": 
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which is coupled to Friedmann's equation, 
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where the final expression is satisfied since we assume the 
Universe is radiation-dominated throughout the phase 
transition. Our choice of £ < 1 allows for this assumption 
to be valid. 

We consider the case where £ = 0.25 and n bubbles 
nucleate simultaneously per Hubble volume, H~ 3 ; our 
simulations will involve a smaller volume, V = H~ 3 /8 
with N = n/8 bubbles. At the time of nucleation, each 
bubble has the field profile of the Coleman bounce [lj| , 



(x, t*) = 0o tanh 




r*) 



where 



(e^oVA)" 



(8) 
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is the initial bubble radius. Although the tunneling rate, 
0, is set by the parameters of the model (A, 4>q and e) 
[16| , we consider a toy model in which the tunneling rate 
is a free parameter, or rather, that the number of bubbles 
per Hubble volume, TV, is a free parameter. 



NUMERICAL SIMULATIONS 

The classical production of gravitational radiation 
from lattice simulations of scalar fields has received a 

TJ-El, particu- 



lot of attention over the last decade 
larly in studying the production of gravitational radiation 
from preheating. The majority of this analysis has em- 
ployed LatticeEasy [26j to simulate scalar fields in an 
expanding background, although different authors have 
chosen different methods of calculating the power in grav- 
ity waves. We continue to use a modified version of Lat- 
ticeEasy for our scalar evolution. 

In general, we will set the energy scale of the simula- 
tions, fx, to 10~ 4 m p i, 10 _5 m p i, and 10 _6 m p i. The other 
parameters in the potential are constrained by £ and by 



where n = 8N is the number of bubbles per Hubble 
volume. After initializing the bubbles of true vacuum, 
we evolve the field in a radiation-dominated expanding 
background until a final time tf ~ H~ 1 /2, taking 10 4 
time steps on a 512 3 grid. Over the course of the sim- 
ulation, the size of the universe increases by a factor of 
a = 1.5. Figs. Q] and [2] show the time evolution of the 
field and gravitational wave spectrum for a simulation 
with p* = (10" 4 m p i) 4 and N — 5 bubbles (n = 40). The 
fields then evolve according to Eq.[S] until t sa 2.5/?~ 1 . 

We use the algorithm of [lj| 21 1 to calculate the 
power spectrum of gravitational radiation produced in 
this phase transition. The classical metric perturbations, 
hij, written in synchronous gauge, are 



ds 2 



dt 2 - a 2, {t) [Sij + hij] dx l dx J . 



(11) 



When we impose the transverse-traceless conditions, 
h TT \ = and h TTl \j = 0, the metric perturbations 
are subject to the equation of motion, 
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The source of this wave equation is the transverse- 
traceless projection of the anisotropic stress tensor, 



Si 



Vij 



T. 



(13) 



We evolve the metric perturbations alongside the scalar 
fields; this allows us to track the power in gravitational 
radiation at any point during the simulation and pin- 
point from where that radiation emerges. The stress- 
energy associated with metric perturbations is [28j, 
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where the brackets denote a spatial average over several 
wavelengths. The 00 component is the energy density, 



Pgw 
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FIG. 1: Four time-slices of a first order phase transition when the energy density of the Universe was p ~ (10~ m p i) . The 
first slice shows the nucleation of five bubbles at r = 0, followed by a slice taken as the bubbles initially collide (r « 0.5/3 -1 ), 
a slice at the end of the phase transition (r ~ and finally a slice when r ~ 2f3~ 1 . Contours are drawn at <f> = — 0.83<^o 

(gold) and <j> = 0.83<^>o (red) to guide the eye. 



which can be converted to momentum space by use of 
Parseval's theorem (see |2lj ) 



™}l 1 
32 



(16) 



where 1/ is the comoving volume over which the spatial 
average is being performed. We then write 



dp. 



g\v 
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m^k 3 1 
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This corres pon ds to a present-day amplitude and fre- 
quency by [2l|, [2J] , 



1/3 
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FIG. 2: Gravitational wave spectra from a first-order phase 
transition when the energy density of the Universe was p w 
(10 _4 m p i) 4 : the spectrum immediately after the nucleation 
of five bubbles (red, solid), at r ~ 0.5/3 - (blue, dotted), 
when r w /3 _1 (green, dashed), and at t w 2f3~ 1 (black, dot- 
dashed). The bump at high frequencies is a numerical artifact. 



/ = 6 X 10 1 



:HZ 



(19) 



where the subscript indicates quantities defined today 
and e can denote any time during (or at the end of) the 
simulation. We also keep the convention that h absorbs 
the uncertainty in the present value of the Hubble pa- 
rameter. The quantity f2 ra( j,o is the current fraction of 
the energy density in the form of radiation, and ptot,c is 
the total energy density at t e . The ratio go/g* compares 
the number of degrees of freedom today to the num- 
ber of degrees of freedom at matter/radiation equality. 
We approximate go/g* = 1/100. Fig. [5] shows four time 
slices (corresponding to the same times as in Fig. Q} of 
the gravitational radiation produced in a simulation with 
p* = (10- 4 m pl ) 4 and N — 5 bubbles (n = 40). 



RESULTS 

As the bubbles expand and begin to collide, individual 
collisions are seen to contribute to the power spectrum. 



These are no longer apparent by t w /? _1 /2, after which 
time h 2 il gw rises more steadily, as shown in Fig. [3] We 
begin by considering the time interval < r < . 



Collisions: < r < /3 1 

During this stage, the peak frequency of the gravita- 
tional radiation should correspond to the scale a e j3~ l — 
L* p}, the mean distance between bubbles on the ini- 
tial slice. This corresponds to a physical wavenumbcr 
fephys = ZkL^ 1 s» 11.34rt 1//3 ii*. The associated frequency 
observed at the present day is given by Eq. [191 where H e 
is the Hubble constant at the time when we calculate the 
spectrum; H e = 4jt ~ H Q , where a e is the scale factor 
at the time when we take the spectrum. Putting this 
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FIG. 3: The maximum intensity of the gravitational wave 
spectrum for a /i = 10 -4 m p \ simulation initialized with 40 
(red, solid), 32 (blue, dotted), 24 (green, dashed), or 16 
(black, dot-dashed) bubbles per Hubble volume, r < /3 _1 . 



together, we expect the peak frequency to occur at 



/peak ~ 6.8 x 10 u 7i 1/3 ,/— Hz 1.16 x 10 12 M7i 1/3 Hz 

(20) 

This corresponds to a peak around / p0 ak ~ 10 8 Hz when 
/i = m p i and n ~ 10. These numbers are consistent 
with the location of a low-frequency peak visible in the 
spectrum at early times; see Fig. 2J 
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FIG. 4: The present-day gravitational wave spectrum pro- 
duced by time r — This is for a simulation with 
H = 10~ 4 m p i and n — 16 bubbles per Hubble volume (red, 
solid), n = 24 (blue, dotted), n — 32 (green, dashed), and 
n = 40 (black, dot-dashed). The bump at high frequencies is 
a numerical artifact. 

The frequency at which the maximum intensity of 
gravity waves is expected to be found is also given in 
[7f as 



peak 



5.2 x lCT 8 Hz 



I 



,(—) 

lGeVy V100/ 



1/6 



(21) 



where g„ is the number of ultra-relativistic degrees of 
freedom at the time of the transition, j3 is, again, the 



nucleation rate, and T* = /u is the energy density when 
the phase transition occurs. Taking g* = 400, 



/p Cak = l.lSn 1 / 3 ^ x 10 12 , 



(22) 



which predicts the spectrum to peak at frequencies three 
or four times \x x 10 12 , depending on the number of bub- 
bles. Fig. [5] confirms that for n = 10~ 4 , peak frequencies 
occur around this value. 

A more recent calculation of the peak frequency 13 1 
predicts 



/ = 16.5 x 10" 



lOOGeVy UOO/ 

(23) 

where we use the same values as before and approximate 
f*/P by 

U _ 0.62 
P ~ 1.8 -Q.lv - 



= 0.22963 



(24) 



when the bubble wall velocity v « 1 (the scalar bub- 
bles in our simulation accelerate quickly to v « 1). This 
prediction of the frequency simplifies to 



/peak = 8.60n 1/3 Ai x 10 



li 



(25) 



which also predicts a few times 10 8 , also consistent with 

Fig. m 

In addition to calculating the peak frequency, also 
gives the fraction of critical density found in gravity 
waves as 



=1.1 x 10 



0.24 



p 

100 

g* 



1/3 



(26) 



where the efficiency factor k measures energy lost to the 
motion of a coupled fluid. We set n = 1, as our simu- 
lations do not include a fluid. The fraction H*/j3 is the 
mean initial bubble separation as a fraction of the Hubble 
distance, v is the velocity of the bubble walls, and g* is, 
again, the number of ultra-relativistic degrees of freedom. 
In our simulations we set a — 1/3. Assuming v = 1 and 
<7* = 400, the expected intensity of the spectrum reduces 
to 



n gw h 2 = 1.42n _2/3 x 10~ 9 . 



(27) 



This estimate varies between 90 and 1250 times greater 
than simulation results at r = P~~ , becoming more con- 
sistent in the large n limit. The more recent prediction 
fl3l | of intensity as 



n gw h 2 =1.67 x 10" 



O.llir 4 
0.42 + 1 

100 

.9* 



a* 
T 



l + a 
=2.07n _2/3 x 10~ 9 



1/3 



(28) 
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falls slightly farther from our results. There are two main 
factors that describe the discrepancies: (1) we have im- 
plemented an expanding background in which Hubble 
friction depletes the energy of the source by a small frac- 
tion and (2) our models realistically thicken the bubble 
walls. "Thick" walls prolong collisions and dilute the gra- 
dient terms that source strong gravitational waves. At 
the same time the inclusion of these effects strengthens 
the validity of the current model. 



Coalescence: /? 1 < r < 2/3 1 : 

Although most of the volume of the simulations is in 
the true minimum by r = there is still a lot of ki- 
netic, gradient and even potential energy in the fields. 
The phase transition is, more or less, complete, but 
the production of gravitational radiation has not ceased. 
Fig. [5] shows how the peak amplitude rises from r = 
until t = 2.5/?- 1 . 




FIG. 5: The maximum intensity of the gravitational wave 
spectrum for a fi — 1CP 4 m p i simulation initialized with 40 
(red, solid), 32 (blue, dotted), 24 (green, dashed), or 16 
(black, dot-dashed) bubbles per Hubble volume, /3 _1 < r < 

2.5/r 1 . 



We expect that a period of turbulence can increase the 
magnitude of the spectrum by several orders of magni- 
tude, e.g. [7-12|. One of these estimates Q say that 
the intensity of gravitational radiation after turbulence 
is predicted to be 



Qgwh = 10 



100 



1/3 



which reduces to 



n gw h 2 = 2.55n" 2/3 x 10~ 7 . 



(29) 



(30) 



Thus, the intensity of the gravitational wave spectrum is 
expected to be around order 10~ 8 if we take Vq « 1. Wc 
can do slightly better if we try to assign a sound speed, 
Vq, by estimating the speed of perturbations in the true 



vacuum. Near <p 
m 2 ff = X(/)q, and 



the effective mass of the field, 



v = 



~dk 



i + \<}>l/k 2 



(31) 



This ranges between 10 -2 for low frequency modes, 
fciow ~ H* s» VA0o/5OO, and almost 1 for higher fre- 
quency modes, k ~ -s/A^o- The final amplitudes that we 
present here, see Fig. [SJ are some three orders of mag- 
nitude lower, but this is understandable as there is no 
turbulence per se in our simulation. 

There is, however, significant post-collisionary ampli- 
fication of the gravitational wave spectrum. This period 
of coalescence after r = of the simulation, ampli- 
fies the spectrum by more than an order of magnitude, 
depositing energy in higher frequency modes. 

The frequency of the turbulence peak is Q 



/ peak ~ 2.6 X IO^HzUoW" 1 



1/6 



H,J VlGeVy V100/ 
= 5.91 x 10 1 V /: V x 10 11 

(32) 

This predicts the peak to shift downward during turbu- 
lence, while we see a higher-frequency peak at the end 
of the simulation, with frequency / ~ O(fi) x 10 13 . This 
peak comes from the amplification of higher frequency 
modes during the coalescence period. Fig. [5] shows the 
integrated gravitational wave spectrum for u. — 10~ 4 af- 
ter this period. In this plot the peak at higher frequencies 
is quite apparent at r = 2.5/3" 1 . 
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FIG. 6: The present-day gravitational wave spectrum at t = 
2fi~ 1 , nearing the end of the coalescence phase. This is for a 
simulation with p» = (10~ 4 m p i) 4 and 16 bubbles per Hubble 
volume H* 3 (red, solid), 24 bubbles per Hubble volume (blue, 
dotted), 32 bubbles per Hubble volume (green, dashed), 40 
bubbles per Hubble volume (black, dot-dashed). The bump 
at high frequencies is a numerical artifact. 

It should be noted that Fig. [5] does not explicitly show 
the k 3 low-frequency tail of the gravitational wave spec- 
trum. This is due to the lack of resolution at the relevant 
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scales; the longest wavelength that we can resolve is well 
within the horizon at the time of the phase transition. 
The leftmost few points in all of our spectra are aver- 
aged only over a small number of modes and should not 
be used to extrapolate to very low frequencies. 

Lastly, it's important to check how the spectrum varies 
with energy scale. In Fig. [7] we vary the energy scale, ^, 
between 10 _6 rn p i and 10~ 4 m p i. The amplitude of the 
gravity wave spectrum should be independent of scale 
of the simulation, fi (also T* in 0, among others), 
and should only depend on /3 _1 and dynamical factors. 
Indeed, we recover the scale- independent behavior for the 
three orders of magnitude that we test. We anticipate 

1 x 10~ 10 i ■ — ■ — ■ — ■ — ■ — . . . . . 



5xi<r" - 




1 x LIT 13 1 1 1 1 1 1 

10 7 10 s 10" 10'" 

f (Hz) 

FIG. 7: The present-day gravitational wave spectra at r = 
2.5/3' 1 for cases with 32 bubbles per Hubble volume at three 
energy scales: fj, = 10 -6 m p i (green, dashed), 10 -5 m p i (blue, 
dashed) and 10 -4 m p i (red, solid). The bump at high frequen- 
cies is a numerical artifact. 

that our simulations would continue to produce similar 
spectra even at much lower energy scales (modified only 
slightly when the number of ultra-relativistic degrees of 
freedom, g e , decreases). 



settles into the true vacuum. During the first stage, we 
precisely reproduce the location of the peak of gravita- 
tional radiation from previous estimates. The amplitude 
at this time is lower than expected from these estimates, 
due to the inclusion of friction and by realistically "thick- 
ening" the walls. 

Primarily, though, we have discovered that a coales- 
cence phase following the phase transition amplifies the 
gravitational wave signal for a first-order phase transition 
by about an order of magnitude and increases the peak 
frequency by about a decade. This is most likely a conse- 
quence of the persistence of energy in domain walls, even 
after regions have collided, e.g. [29|], residual anisotropic 
stress-energy produced as the Universe relaxes to a ther- 
mal state. This compensates for the lack of power in 
gravity waves at t s» /3 _1 . 

The electroweak phase transition will occur at a much 
lower frequency than those presented here. If we esti- 
mate this energy scale as ~ 200 GeV, then we expect 
the peak frequency of gravitational radiation from the 
phase transition to occur at a few times 10 -5 Hz, assum- 
ing P/H* « 5 as we have here. Coalescence will both 
amplify the signal and raise the peak frequency about an 
order of magnitude. We believe that this is an important 
effect to be considered in the next generation of detector 
experiments. 

We intend to follow up on this model by adding dynam- 
ical fluids to our 3-dimensional simulations. The evolu- 
tion of fluids alongside our scalar fields will allow us to 
confirm the parametric dependence of the gravitational 
wave signal on the terminal velocity of the bubble walls, 
check for the amplification of the signal due to the exis- 
tence of turbulence and confirm that a coalescence phase 
exists in the presence of a viscous fluid. 
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DISCUSSION 

Gravitational radiation should be the most obvious 
relics of first-order phase transitions that may have ex- 
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the predicted scalings of the results both in frequency 
and amplitude. We had identified the two relevant stages 
of the process: (I) the stage during which the bubbles 
collide and (2) a coalescence phase during with the field 
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